#ifndef MODEL_
#define MODEL_

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <cmath>
#include "ANN/ANN.h"
#include <omp.h>

#include "DGMatrix.h"
#include "DGIntMatrix.h"
#include "mesh.h"
#include "saveGmsh.h"
#include "initialSol.h"
#include "functionSpace.h"
#include "flux.h"
#include "muparser/muParser.h"

using namespace std;

//enum PType    {EQP,ENP,FQP,FNP,BQP,BNP};
//------------------------------------------------------------------------------
class Model{

private:
  
  int  order;
  bool graphLap;
  bool lumpedM;
  
private:
    
  SolContainer UEQP;                                                            // [EType] (nbEQP,nbE)
  SolContainer UFQP;                                                            // [L&R]   (nbFQP,nbF)
  SolContainer UBQP;                                                            // [L&R]   (nbFQP,nbB)

  SolContainer UFNP;                                                            // [L&R]   (nbFNP,nbF)

  vector<SolContainer> FUEQP;                                                   // [EType] [dim] (nbEQP,nbE)
  vector<SolContainer> FUFQP;                                                   // [L&R]   [dim] (nbFQP,nbF)
  vector<SolContainer> FUBQP;                                                   // [L&R]   [dim] (nbFQP,nbB)
  
  vector<SolContainer> FUFNP;                                                   // [L&R]   [dim] (nbFNP,nbF)
  
  vector<SolContainer> dFUEQP;                                                  // [EType] [dim] (nbEQP,nbE)
  vector<SolContainer> dFUFQP;                                                  // [L&R]   [dim] (nbFQP,nbF)
  vector<SolContainer> dFUBQP;                                                  // [L&R]   [dim] (nbFQP,nbB)

private:
  
  SolContainer bE;                                                              // [eType] (3,3) --> elmt bilinear form
  DGMatrix     bF;                                                              //         (4,4) --> face bilinear form
  
  vector<DGVector> nuE;                                                         // [eType] (eID)
  DGVector         nuF;                                                         // (fID)
  
  void setBilinearForm                ();
  void computeViscosity               (const SolContainer*);
  void computeArtViscousVolumeTerm    (const SolContainer*,SolContainer*);
  void computeArtViscousInterfaceTerm (const SolContainer*,SolContainer*); 
  
private:
  
  void computeFUQP           ();                                                // compute flux at quadrature points (elements and faces)

  void computeVolumeTerm     (SolContainer*);
  void computeInterfaceTerm  (SolContainer*);
  void computeBoundaryTerm   (SolContainer*);

  void multiplyMass          (SolContainer*);

public:
  
  Flux*           flux;
  Mesh*           mesh;
  FsContainer*    fs;
  SaveOutputGmsh* output;

public:
  
  Model(Mesh*,FsContainer*,SaveOutputGmsh*,ClassParameter&);
  ~Model(){};
  
  void computeResidual      (const SolContainer*,SolContainer*);
  void computeUQP           (const SolContainer*);                              // compute solution at quadrature points (elements and faces)
  void computeMaxPrinciple  (const SolContainer*,DGMatrix*);  

  const SolContainer* getUEQP() {return &UEQP;}
  
  double getCFL(double);

  
//---ici---
/*  
  
private:
  
  DGMatrix* bE;
  DGMatrix* bF;
  
  map<int,vector<double> >  nuE;
  vector<double>            nuF;
  vector<double>            nuBC;
  
  void setBilinearForm();
  
private:
  
  void getLinearConvectiveFlux        ();
  void computeBoundaryTerm            (SolContainer*);
  void computeViscosity               (SolContainer*);
  void computeArtViscousVolumeTerm    (SolContainer*);
  void computeArtViscousInterfaceTerm (SolContainer*);  

public:
  
  void computeResidual      (SolContainer*,SolContainer*);
  void collocateSolution    (SolContainer*);
  
  double getCFL(double);
*/  
};
#endif
